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Abstract 

A population of firing neurons is expected to carry not only mean firing rate 
but also its fluctuation and synchrony among neurons. In order to examine 
this possibility, we have studied responses of neuronal ensembles to three kinds 
of inputs: mean-, fluctuation- and synchrony-driven inputs. The generalized 
rate-code model including additive and multiplicative noise (H. Hasegawa, 
Phys. Rev. E 75, 051904 (2007)) has been studied by direct simulations (DSs) 
and the augmented moment method (AMM) in which equations of motion for 
mean firing rate, fluctuation and synchrony are derived. Results calculated 
by the AMM are in good agreement with those by DSs. The independent 
component analysis (ICA) of our results has shown that mean firing rate, 
fluctuation (or variability) and synchrony may carry independent information 
in the population rate-code model. The input-output relation of mean firing 
rates is shown to have higher sensitivity for larger multiplicative noise, as 
recently observed in prefrontal cortex. A comparison is made between results 
obtained by the integrate-and-fire (IF) model and our rate-code model. The 
relevance of our results to experimentally obtained data is also discussed. 
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1 Introduction 



One of the most important and difficult problems in neuroscience is to understand 
how neurons communicate in a brain. There has been a long-standing controversy 
between the temporal- and rate-code hypotheses in which information is assumed to 
be encoded in firing timings and rates, respectively [1]- [3]. A recent success in brain- 
machine interface (BMI) [lj[5], however, suggests that the population rate code is 
employed in sensory and motor neurons, though it is still controversial which of rate, 
temporal or other codes is adopted in higher-level cortical neurons. 

In recent years, much attention has been paid to a study on effects of mean fir- 
ing rate, its fluctuation and synchrony (or spatial correlation) of input signals (for a 
review on rate and synchrony, see Ref . j6] ) . The precise role of synchrony in informa- 
tion transmission and the relation among the firing rate, fluctuation and synchrony 
are not clear at the moment [6]- [12]. The firing rate and synchrony are reported 
to be simultaneously modulated by different signals. For example, in motor tasks 
of monkey, firing rate and synchrony are considered to encode behavioral events 
and cognitive events, respectively [TJ. During visual tasks, rate and synchrony are 
suggested to encode task-related signals and expectation, respectively [8] . A change 
in synchrony may amplify behaviorally relevant signals in V4 of monkey [9]. An in- 
crease in synchrony of input signals is expected to yield an increase in output firing 
rate. The synchrony of neurons in extrastriate visual cortex is, however, reported to 
be modulated by selective attention even when there is only small change in firing 
rate [12]. Rate-independent modulations in synchrony are linked to expectation, 
attention and livalry [6] . Fluctuations of input signals have been reported to modify 
the / — / relation between an applied dc current / and autonomous firing frequency 
/ although its sensitivity to input fluctuation seems to depend on a kind of neurons 
[13j-[15j. The / — I curve of prefrontal cortex (PFC) retains the increased sensi- 
tivity to input fluctuations at large J, while that of somatosensory cortex (SSC) is 
insensitive to input fluctuation though its linearity is increased at small / [T5] . 

This kind of problems discussed above have been extensively studied by using 
spiking neuron models such as the Hodgkin-Huxley model [H] and integrate-and-fire 
(IF) model with diffusion and mean-field approximations [UJ-[33]. The purpose of 
the present paper is to examine the same problem by using the rate-code model, 
which is an alternative theoretical model to the spiking model. In a previous paper 
[51] (referred to as I hereafter), we proposed the generalized rate-code model for 
coupled neuron ensembles with finite populations, which is subjected to additive 
and multiplicative noises. It seems natural to include multiplicative noise beside 
additive noise in our rate-code model because the noise intensity is expected to 
generally depend on the state of neurons. Actually effects of multiplicative noise 
in the spiking neuron model are extensively examined by using conductance-based 
inputs which yield multiplicative noise [32]" [HZ]- Our calculation in I has shown that 
the introduced multiplicative noise leads to the non-Gaussian stationary distribu- 
tion of firing rate, yielding interspike-interval (ISI) distributions such as the gamma, 
inverse-Gaussian and log-normal distributions, which have been experimentally ob- 
served. We have discussed the dynamical properties of neuronal ensemble, by using 
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the augmented moment method (AMM) which was developed for a study of stochas- 
tic systems with finite populations [38j • In the AMM, we pay our attention to global 
properties of neuronal ensembles, taking account of mean, and fluctuations of lo- 
cal and global variables. The AMM has the same purpose to effectively study the 
properties of neuronal ensembles as approaches based on the population-code hy- 
pothesis [3H]-[I3]- The AMM has been nicely applied to various subjects of neuronal 
ensembles [44j [45] and complex networks [46J. 

We have assumed in I that input signals are the same for all neurons in the 
ensemble. In the present study, input signals are allowed to fluctuate and to be 
spatially correlated. We will derive equations of motion for mean firing rate, its 
fluctuation and synchrony with the use of the AMM, in order to investigate the 
response to mean-, fluctuation- and synchrony- driven inputs. This study clarifies, 
to some extent, their respective roles in information transmission. 

The paper is organized as follows. In Sec. 2, we discuss an application of the 
AMM to the generalized rate model, studying the input-output relations of mean, 
fluctuation and synchrony. In Sec. 3, stationary and dynamical properties are 
discussed with numerical model calculations. By using the independent component 
analysis (ICA) [47J, we investigate a separation of signals when two or three kinds 
of inputs are simultaneously applied. Discussions are presented in Sec. 4, where 
the results obtained by spiking IF model are compared with those by our rate-code 
model. The final Sec. 5 is devoted to conclusion. 



2 Formulation 

2.1 Adopted model 

For a study of the properties of a neuron ensemble containing finite N neurons, we 
have adopted the generalized rate-code model [34"1 14"5] in which a neuron is regarded 
as a transducer from input to output signals, both of which are expressed in terms 
of firing rates. The dynamics of the firing rate r^(t) (> 0) of a neuron i (i — 1 to 
N) is given by 



dt 

with 



dv ■ 

= F{r i ) + H{u l ) + G{ ri ) m {t)+ &(t), (1) 



«*(*) = (Vl E rAV+m, (2) 

1 z J m) 



ii 

H(u) = -_=e(n). (3) 
V w + 1 

Here F(r) and G(r) are arbitrary functions of r, Z (= N — 1) denotes the coordina- 
tion number, 7j(t) an input signal from external sources, w the coupling strength, 
and Q(u) the Heaviside function: 0(it) = 1 for u > and Q(u) — otherwise. Addi- 
tive and multiplicative noises are included by £j(t) and rji{t), respectively, expressing 
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zero-mean Gaussian white noise with correlations given by 



(Vi{t)ViM) = a%5(t-t'), (4) 
= P%At-t'), (5) 
(fkWtiffl) = 0, (6) 

where the bracket (■) denotes the average over the distribution of p({ri},t) [Eq. 
( IA3I) ]. and a and (3 denote the magnitudes of multiplicative and additive noise, 
respectively. The gain function H(u) in Eq. expresses the response of a rate 
output (r) to a rate input (u). It has been shown that, when spike inputs with 
mean firing rate are applied to a Hodgkin-Huxley (HH) neuron, mean firing rate 
r a of output signals is r Q ~ for ~ 60 [Hz], above which r a shows the saturation 
behavior [IEJ HH]. The nonlinear, saturating behavior in H(u) arises from the fact 
that a neuron cannot fire with the rate of r > l/r r (= r max ) where r r denotes the 
refractory period. The function H(u) has the rectifying property because the firing 
rate is positive, which is expressed by the Heaviside function in Eq. (j3J). Although 
our results to be present in the following are valid for any choice of H(x), we have 
adopted, in this study, a simple analytic expression given by Eq. (jHJ), where the rate 
is normalized by r max . 

With the use of the diffusion-type approximation, a spatially correlated input 
signal Ii(t) in Eq. (j2J) is assumed to be given by 

ii(t) = mW + siiit), (7) 

with 

(sim = o, (8) 

(SmSIjit')) = [6 ijlI (t) + (l-6 lJ )C I (t)]5(t-t'), (9) 
where variance (7/) and covariance ((j) are given by 

u(t) = ^£<«<(f) a ), (10) 

i 

We will discuss responses of the neuronal ensemble described by Eqs. ([I])-© to the 
spatially correlated input given by Eq. using both direct simulations (DSs) and 
the AMM ElEHl. 



2.2 Augmented moment method 

In the AMM [38J , we define the three quantities of /1, 7 and p given by 

Mt) = W)> = -^ !>*(*)>. ( 12 ) 
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7(*) = ^EW*)-M*)] a >, (13) 

= 4eE(W^)-M*)][^)-^)]>> (1 4 ) 

« 3 

where R = (l/N)J2i r i, (i>{t) expresses the mean, and j(t) and p(t) denote the 
averaged, auto and mutual correlations: related discussion will be given later [Eqs. 

flUD-n]. 

By using the Fokker-Planck equation (FPE), we get equations of motion for (rj) 
and (r^j) (for details see appendix A). Expanding r t in Eqs. (|A40 and (IA50 around 
the average value of p as 

r i = p J r 8r h (15) 

and retaining up to the order of (SriSrj), we get equations of motion for p, 7 and 
p. AMM equations in the Stratonovich representation are given by (the argument t 
being suppressed) 



dp 
~dt 



fo + h + M 



'a 2 ' 



^ = 2/ l7 + ^ (Ap - 7 ) + 7J + 2^ + 2< M2 )a 2 7 + a 2 ^ 2 + (3 2 , (17) 
^ = 2f lP + 2h lW p+ + Z(j) +2(g 2 1 +2g g 2 )a 2 p+ ^(a 2 g 2 + (3 2 ),(18) 

where £ = (l/£!)(<^F(p)/«9x £ ), g t = (l/£\)(d e G(p)/dx e ), h t = {l/£\)(&H(u)/du?) 
and u = wp + pj. Original A-dimensional stochastic DEs given by Eqs. (til)-© 
are transformed to the three-dimensional deterministic DEs given by Eqs. (1TB]) - 
( TTSI) . For 7/ = (j = 0, equations of motion given by Eqs. (fTBjl - tfTS]) reduce to those 
obtained in our previous study [31]. From p, 7 and p obtained by Eqs. (fT6]) - (fl8]) . we 
may calculate important quantities of synchrony and variability. Then, they may be 
expressed in physically more transparent forms, as will be discussed in the following 

[Eqs. Q2BD-PBD]. 



2.3 Rate synchronization and variability 
2.3.1 Synchronization ratio 

The synchronization is conventionally discussed for firing timings (temporal synchro- 
nization) or phase (phase synchronization) . We discuss, in this paper, the synchro- 
nization for firing rate (rate synchronization) . In order to quantitatively discuss the 
synchronization, we first consider the quantity P(t) given by 

Pit) = ^ E < Ht) - r,(t)} 2 >= %(*) - Pit)]- (19) 

ij 
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When all neurons are firing with the same rate (the completely synchronous state), 
we get Ti(t) = R(t) for all i, and then P(t) = in Eq. (TT9j) . On the contrary, we get 
P(t) = 2(1 - l/iV)7(t) = P (t) in the asynchronous state where p = 7/JV [31 [38]. 
We may define the normalized ratio for the synchrony of firing rates given by 



where 



s(t) = 1 - *® - ( MlMLzl ) - m m 
S{t) - 1 Po(t)-{ N-i )- l{ ty (20) 
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VZ EEWW))- ( 22 ) 

S'(i) is and 1 for completely asynchronous (P = Pq) and synchronous states 
(P = 0), respectively. 

A similar analysis yields that the synchrony ratio of input signals may be ex- 
pressed by 

Sl - N-^.asm ~v (23) 

where 7/ and £i are given by Eqs. (ITUj) and ffTTl) . respectively. 
2.3.2 Variability 

The variability in the ISI is usually defined by by Cy = \fi^ / p T where /i T and j T 
stand for mean and variance of ISI, respectively. We here define the rate variability 
given by 

Cy[t > ~ — -p. — -755- (24) 

Similarly, the variability in input rate signals is given by 



Quantities calculated in the rate-code hypothesis are compared to those obtained 
in the temporal-code one in the Appendix C. We show in Eq. flOlTj) that if the 
temporal firing rate is given by r^t) = 1/Tj(i), we get Cy(t) — Cy(t), where T^t) 
denotes the interspike interval (ISI) of firing times in a neuron %. 



6 



2.4 AMM equations for fi, 7 and S 

It is noted that the variability and synchrony are related to the second-order statis- 
tics of local (7) and global fluctuations (p), respectively, of firing rates. Employing 
the relations given by Eqs. (1201) and (I23p . we may transform equations of motion 
for p, 7 and p given by Eqs. (fTB]) - (fI5|) to those for p, 7 and S given by 

dp (a 2 \ 

= fo + h + / 2 7 + ( y J [9o9i + 3(#i# 2 + Qogsh], (26) 

^ = 2fa + 2h 1 w 1 S + 2(g 2 1 +2g g 2 )a 2 1 + ll + a 2 g 2 +P 2 J (27) 

f = ~(0 + «% 2 + ^ + (^) + (^)(l + ^)(l- 5 ) ; (28) 

where f e = (l/£\)(d e F(p) /Ox 1 ), g t = (l/£\)(d i G(p)/dx e ), h e = (l/£\)(d e H(u)/du e ) 
with u = wp+pj. Equations (T26]) - (T28|) express a response of (p, 7, S) to a given input 
of (pi,ji, Sj). Input signals in which information is encoded in pi, 77 and Si are 
hereafter referred to as mean-driven, fluctuation-driven and synchrony-driven inputs, 
respectively. Equations (!26|) - (|28|) show that 7/ plays a similar role to independent 
noise while ( T (= 77-SV) to correlated noise. 
When we adopt F(r) and G{r) given by 

F(r) = -Ar a , (29) 
G(r) = r\ (a,6>0) (30) 

Eqs. fl26D-([28D become 



g?7 
~dt 



+ (j) + b (b - 1)(26 - l)/i ib - 3 7], (31) 

= -2Aa / u a " 1 7 + 2/iiW7S + 26(26- l)aV fe " 2 7 + 7/ + «V b + /^ 2 , (32) 



f = -W + «V fc + /5 2 )5+(^) 

+ (^) (1 + ^(1-5), (33) 

where A expresses the relaxation rate. We note that for (i) a = or 1 and (ii) 6 = 0, 
1/2 or 1, a motion of p is decoupled from the rest of variables because the a(a — 1) 
and b(b — l)(2b— 1) in the third and fourth terms of Eq. (13T|) vanish. Equation (j33l) 
shows that a motion of S is ostensibly independent of the index a of F(r) although 
it depends on a through 7. 
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2.5 Relevance to experimental data 

Data obtained by neuronal experiments have been analyzed by using various meth- 
ods [6] [50]- [52]. It is worthwhile to discuss the relevance of experimentally observed 
data to j(t), p(t), ((t) and S(t) which are introduced in our study. Let x^(t) be a 
binned data defined by 

x1(t) = 1, if 4 e [t- At/2, t + At/2) 

= 0, otherwise (34) 

where t k in stands for the nth firing time of neuron i of a given kth trial and At the 
width of time window. The firing rate of the neuron i averaged over M trials is 
given by 

/ I \ M 

r <MMAi)l>- (t) - (35) 
One conventionally calculates the correlation averaged over time given by 

/oo 
T i3 {t,t + r)dt, (36) 
-oo 

with 

r-(M') = MOMO-MOl, (37) 

where «(t) denotes the average of firing rates given by 

1 N 

/#) = i ( 38 ) 
iy i=i 

However, Cy(r) is not so useful for a study of the response to time- dependent 
inputs as in our case. The equal-time, auto and mutual correlations averaged over 
an ensemble are given by 

A (t) = ^E r *(M), (39) 

i 

M(t) = ^EE r «(M). (40) 

Comparing Eqs. flU and ([22]) with Eqs. ([39]) and gO]), we get 

l(t) = A(t), (41) 

C(t) = M(t), (42) 

which yields 

p(t) = L[A(t) + (N-l)M(t)}, (43) 

m - m 

Equations fl4"TT) - (l44"I) show that j(t) and ((t) are equal-time, auto and mutual corre- 
lations, respectively, and that S(t) is nothing but the mutual correlation normalized 
by the auto correlation. 
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3 Model Calculations 



3.1 Stationary properties 

In order to get an insight to the present method, we will show some model calcula- 
tions. When we consider a special case of a = b = 1.0 in Eqs. ( 1291 and ( 1301 : 



F(r) = -Ar, (45) 
G(r) = r, (46) 



Eqs. fl3Tl) - fl33l) are expressed by 



i = + < 47 > 

-2 = -2A7 + 2h lW ^S + 7/ + 2a 2 7 + a V + P 2 , (48) 
at 

f = -^(2c 7 , + aV + /3 2 ) + ^+(^)(l + ^)(l-S), (49) 
The stationary solution of Eqs. (I47l) -(l49l) is given by 

" = (A-aV2) ' (50) 

( 7/ + a 2 /i 2 + /? 2 + 2h lW Np) 

7 2(A - a 2 + h lW /Z) ' 1 J 

( 7/ + aV + /? 2 )[Z(A-a 2 ) -h lW (Z- 1)] + h^ I S I ' 1 J 

" ;/ ' S/ forw = (53) 



(7/ + a V + /3 2 



[Z(A -a 2 ) - h lW (Z~l)] 



for 5/ = (54) 



We note in Eq. (!50l that /1 is increased as u/ is increased with an enhancement 
factor of 1/(A — a 2 /2), but « is independent of 7/ and 5/. Local fluctuation 7 is 
increased with increasing input fluctuation (7/) and/or noise (a, (3) as Eq. ( I3T1) 
shows. Equation fl52|) shows that 5 is increased with increasing Sj, as expected. 

The stability condition around the stationary state may be examined from eigen- 
values of the Jacobian matrix of Eqs. (j47p - (|49p . which are given by (for details, see 
appendix B) 



a 2 



Ai = -A + — + hiw, (55) 



A 2 = -2A + 2a 2 -^, (56) 
A 3 = -2A + 2a 2 + 2htw. (57) 
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The first eigenvalue of Ai arises from an equation of motion for p, which is decoupled 
from the rest of variables. The stability condition for p is given by 



hiw < (A - a 2 / 2 )- 



(58) 



In contrast, the stability condition for 7 and p is given by 



Z(\-a 2 ) < htw < (A - a 2 ). 



(59) 



Then for A — a 2 < h±w < A — « 2 /2, 7 and p are unstable but p remains stable. 

The parameters in our model are a, /3, w and N: we hereafter set A = 1.0. Input 
signals are characterized by pi, 7/ and Sj. We will present some numerical calcu- 
lations for mean-, fluctuation- and synchrony- driven inputs applied to an ensemble 
with a = 0.0, j3 — 0.1 and N = 100 otherwise noticed. 

A. Mean-driven inputs 

Figures 1 shows the pi dependences of p and 5 for w = 0.0 (dashed curves) and 
w = 0.5 (solid curves) with 77 = 0.2 and Si = 0.2. We note that for w = 0, p is 
increased with increasing pi after the gain function of H(u), while S is independent 
of pi. In contrast, for finite w, p is much increased than that for w = 0. The chain 
line expressing p = pi crosses the p curve at p = for w = 0.0 and at p = 0.735 
for w = 0.5, 

B. Fluctuation-driven inputs 

Figure 2 shows Cy and S against Cyi (= ^/ji/pi) for w = 0.0 (dashed curves) 
and w = 0.5 (solid curves) with pi = 0.2 and Si = 0.2. For w = 0.0, an increase in 
Cyi yields an increase in Cy and S, while p is independent of CVi, as Eqs. fl50l) -(l52l) 
show. The chain curve shows Cy = Cyf. the region where Cy < Cyi is realized for 
Cyi > 0.51 for w = 0.0 and C VI > 0.22 for w = 0.5. 

C. Synchrony-driven inputs 

Figure 3 shows the Si dependences of p and 5* for w = 0.0 (dashed curves) and 
w = 0.5 (solid curves) with pi = 0.2 and 7/ = 0.2. With increasing Si, S is increased 
as expected. For w — 0, p, is independent of Si, as Eqs. fl50|) shows. For w = 0.5, S 
is much increased compared to that for w = 0.0. 

It is necessary to point out that the pi dependence of p is modified by multi- 
plicative noise (a). An example of the pi dependence of p is plotted in Fig. 4 for 
various a. With increasing a, p shows a steeper increase for larger a because of the 
(A — « 2 /2) factor in Eq. (!50|) . This reminds us the recent experiment of prefrontal 
cortex (PFC) showing that the f — I curve has the increased sensitivity at large 
/ with increasing input fluctuation [15]. This is interpreted as due to a shorten, 
effective refractory period by fluctuation in the calculation using the IF model [15] . 

The dependence of S on Si is plotted in Fig. 5 for various values of N , ensemble 
size (a = 0.5, (3 = 0.2, w = 0.5). It is shown that the synchrony 5* is more increased 
in smaller system. The result for = 100 is nearly the same as that for N = 00. 
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3.2 Dynamical properties 
3.2.1 Pulse inputs 

In order to study the dynamical properties of the neuronal ensemble given by Eqs. 
CO)-©, we have performed direct simulations (DSs) by using the Heun method 
[53| [51] with a time step of 0.0001: DS results are averages of 100 trials otherwise 
noticed. AMM calculations have been performed for Eqs. (j4"Tj) - (|49p by using the 
fourth-order Runge-Kutta method with a time step of 0.01. We consider, as a 
typical example, an ensemble with A = 1.0, a = 0.1, (3 = 0.1, w = 0.5 and N = 100. 
Calculated responses to mean-, fluctuation- and synchrony-driven pulse inputs are 
shown in Figs. 6(a), 6(b) and 6(c), where solid and dashed curves show results of 
the AMM and DS, respectively. 

A. Mean-driven inputs 

First we apply a mean-driven pulse input given by 

fjt^t) = A Q(t - 40)6(60 — t) + A b , (60) 

with A = 0.4, A b = 0.1, 7/(t) = 0.1 and SV(t) = 0.1. Four panels in Fig. 6(a) show 
[A, 7, S and Cy\ chain curves express input signals of 7/(£) and Sj(t). An 

increase in an applied mean-driven input at 40 < t < 60 induces an increase in 
and decreases in j(t) and S(t) which arise from the h\ term in Eqs. f|4*9|) . By an 
applied pulse input, Cy(t) is decreased because of the increased /i. The results of 
the AMM are in fairly good agreement with those of DS. 

B Fluctuation-driven inputs 

Next we apply a fluctuation-driven input: 

7j (t) = B 6(t - 40)9(60 - t) + B b , (61) 

with B = 0.2, B b = 0.05, ///(£) = 0.1 and Si(t) = 0.1, which are plotted by chain 
curves in Fig. 6(b). When the magnitude of ji(t) is increased at 40 < t < 60, j(t) 
and Cy{t) are much increased, while there is no changes in fj,(t). S(t) is modified 
only at t ~ 40 and t ~ 60, where the input pulse is on and off. 

c. Synchrony-driven inputs 

We apply a synchrony-driven input: 

S^t) =Ce(t- 40)0(60 -t) + C b , (62) 

with C = 0.4, Cb = 0.1, Hiif) = 0.1 and ji(t) = 0.1, which are plotted by chain 
curves in Fig. 6(c). An increase in synchrony-driven input at 40 < t < 60 induces 
increases in S(t), j(t) and Cy{b), but no changes in fi(t). This is because fj,(t) is 
decoupled from the rest of variables in Eqs. ( )471) - (l49l) for the case of F(r) = — Ar 
and G(r) = r. 
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3.2.2 Sinusoidal inputs 

We will apply mean-, fluctuation- and synchrony-driven sinusoidal signals to a typ- 
ical ensemble with a = 0.1, j3 = 0.1, w = 0.5 and iV = 100. 

A. Mean-driven inputs 

First we apply a mean-driven sinusoidal input given by 

///(t) = 0.2 [1 - cos(27rf/20)] + 0.1, (63) 

with ji(t) = 0.1 and Si(t) = 0.1. Four panels in Fig. 7(a) show /u(t), 7(t), S(t) and 
Cyit) calculated by the AMM (solid curves) and DS (dashed curves): chain curves 
express input signals of ///(i), Ji{t) and Si(t). An applied signal induces sinusoidal 
output in /i/(t), and small changes in other quantities, as in the case of pulse inputs 
shown in Fig. 6(a). 

B. Fluctuation-driven inputs 

Next we apply fluctuation-driven sinusoidal inputs given by 

7/ (t) = 0.05 [1 - cos(2vrt/20)], (64) 

with Hiit) = 0.1 and Sj{t) = 0.1. Calculated results are shown in Fig. 7(b), 
whose four panels show that an applied input induces sinusoidal changes in j(t) and 
Cy(t), but no changes in /i(t). S(t) shows a peculiar time dependence which may be 
compared to that in the case of pulse inputs shown in the third frame of Fig. 6(b). 

B. Synchrony-driven inputs 

When we apply the synchrony- driven sinusoidal input given by 

Sj(t) = 0.25 [1 - cos(27rt/20)], (65) 

with = 0.1 and ji(t) = 0.1. it leads to a significant change in S(t) and also 

small changes in j(t) and Cy(t), but no changes in fi(t), as in the case of pulse 
inputs shown in Fig. 6(c). 

3.3 Independent component analysis 

It is interesting to estimate multivariate input signals from multiple output signals. 
Such a procedure has been provided in various methods such as Bayesian estimation 
and independent component analysis (ICA) [17]. Here we consider ICA, which was 
originally developed for a linear mixing system, and then it has been extended 
to linear and nonlinear dynamical systems. ICA has revealed many interesting 
applications in various fields such as biological signals and image processing. A 
vector x of output signals is a real function F of a vector s of input sources: 

x=F(s). (66) 

The dimension of s is assumed to be the same or smaller than that of x. If components 
of s are statistically independent and if only one of the source signals is allowed to 
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have a Gaussian distribution, ICA may extract a vector y with a function G given 
by 

y=G(x), (67) 
from which we may estimate the original source as s ~ y [47J . 

A. Coexistence of f^j and Si 

We will discuss the case when mean- and synchrony- driven inputs are simultane- 
ously applied to the neuronal model. We consider the mean-driven sinusoidal input 
and synchrony-driven toothsaw input, given by 

///(t) = 0.1 [1 - cos(27rf/20)] + 0.1, (68) 
Sj(t) = 0.01mod(t,50), (69) 

with 7/(t) = 0.1 where mod(a, b) denotes the mod function expressing the residue of 
a divided by b. Two panels in Fig. 8(a) show /i/(t) and Si(t), and those in Fig. 8(b) 
show ji{t) and S(t) calculated by the AMM. We note a little distortion in S(t) due 
to a cross talk from //(£). Assuming s = (///, S/)t and x = (fj,,S)^, we have made an 
analysis of our result by using ICA. Two panels in Fig. 8(c) show two components 
of y extracted from x = (/x, SY shown in Fig. 8(a) with the use of the fast ICA 
program [55]. Although the program is designed for linear, mixing signals, we have 
employed it for our qualitative discussion. We note that results in Fig. 8(c) fairly 
well reproduce the original, sinusoidal and toothsaw signals shown in Fig. 8(a). 

B. Coexistence of ///, 77 and Si 

Next we consider the case where three kinds of inputs are simultaneously ap- 
plied. They are mean-driven sinusoidal signal, fluctuation-driven toothsaw signal 
and synchrony-driven square pulse signal, given by 

fju(t) = 0.1 [1 - cos(27rf/20)] + 0.1, (70) 
7 7 (t) = 0.002 mod(t, 50), (71) 
Si(t) = 0.5 6(-cos(27rt/120)). (72) 

Three panels in Fig. 9(a) show the input signals of fii(t), 7/(t) and S/(t). Output 
signals of n(t), y(t) and S(t) calculated in the AMM are shown in three panels of 
Fig. 9(b). 7(t) and S(t) are a little distorted by a cross talk. We have made an 
analysis of our result by using ICA, assuming s = (///, 77, SiY and x = (p,y,S)^. 
Three panels of Fig. 9(c) show signals extracted by ICA. Extracted sinusoidal and 
square signals are similar to those of input signals, though the fidelity of a toothsaw 
signal is not satisfactory. This is partly due to the fact that the fast ICA program 
adopted in our analysis is developed for linear mixing models, but not for dynamical 
nonlinear models [55] . 

These ICA analyses show that the mean rate, fluctuation and synchrony may 
independently carry information in our population rate-code model. 
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4 Discussion 



4.1 A comparison with previous studies 

Various attempts have been proposed to obtain the firing-rate model, starting from 
spiking neuron models [56j-[59j. It is difficult to analytically calculate the firing rate 
based on the firing model, except for the IF- type model [l9j-[33]. In the coupled IF 
model, the dynamics of the membrane voltage Vi(t) of the neuron i (= 1 — N) is 
given by 

§ = ~ +m, (73) 

OA, T m 

where r m denotes the relaxation time of the membrane and Ii(t) stands for an input 
to neuron i. When the mean-field and diffusion approximations are adopted, the 
input signal is given by [IS] 

I l (t) = J N (t)+P(t)Ut), (74) 



where (3{t) = J J^fi(t), J is the all-to-all coupling, /i/(t) and jj(t) denote the mean 
and fluctuation, respectively, of input signals, and £j(t) expresses zero-mean Gaus- 
sian white noise with correlations given by (£i(t) £j(t')) = 5ij 5(t — t'). The firing of 
neuron is assumed to occur at t — tf when the voltage v(t) crosses the threshold 
9 from below, and then the voltage is reset to the potential v r : v(tf) = 6, and 
v(tf + 0) — v r . The firing rate r(t) is expressed by 



rif) = ^! 



dp(v, t) 
dv 



(75) 



where the probability distribution p(v, t) is calculated by the FPE with the boundary 
conditions at v = 9 and v = v r and the normalization condition. For the stationary 
state, we get 



I = Tr+ (jL ) f 6 due ^(u-i/x m f/f f u dve -x m (v-i/x m) ye^ (76) 



where / = J/ij, X m =l/r m and r r stands for the refractory period. In order to 
discuss the dynamics of firing rate, it is necessary to solve the FPE to get p(v, t) 
by numerical methods, as was made in Ref. [23]. Equation (|75l) shows that if 
information of input signal is encoded in fluctuation of Ji(t), its transmission is 
instantaneous, as experimentally observed [14], 123]. 

By adopting the IF model, Renart et al. [32] have heuristically derived effective 
equations of motion for the average firing rate \iy and variance of the ISI Oy given 
by 

dfJ-v 1 /™s 

—rr = Vv + Vs, (77) 

at T m 

dvy 2 2 2 

~dT = -^ av + a - (78) 
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where fi s and a s are determined from the stationary solution obtained from the 
FPE of the IF model [as given by Eq. (j7B]) ]. We note that Eqs. ([77]) and JSJ) are 
equivalent to Eqs. (14"7j) and (1481) for A = l/r m and a = 0. Our AMM calculation 
provides us with not only equations of motion for mean and fluctuation but also 
that for synchrony. 

Calculations with the use of the IF model have shown the followings: 

(1) increased input firing rate deceases output variability [32 J, 

(2) increased input firing rate decreases synchrony [20l [22] , 

(3) increased input fluctuation raises firing rate [T5l [26], 

(4) increased input synchrony increases firing rates [H [T2J, [181 E21 El] , and 

(5) increased input synchrony increases variability |21j . 

The items (1), (2) and (5) are consistent with our result shown in Figs. 6(a) and 
6(b). In contrast, items (3) and (4) seem to inconsistent with our result showing 
that fi(t) is independent of ji(t) and Sj(t), as given by Eqs. (1471) . It is noted that 
the model calculation given in the preceding section has been made for the case of 
F(r) = —Xr and G(r) = r, in which a motion of /i(i) is decoupled from those of 
7(t) and S(t). This is not the case in general. We note in Eq. (13T1) that fi{t) is not 
decoupled from y(t) and S(t) except for the case in which (a = or 1) and {b = 0, 
1/2 or 1). For example, in the case of F(r) = —Xr and G(r) = r 2 , equations of 
motion given by Eqs. ( I3~l~l) -(l33l) become 



~dt 

g?7 

~dt 

dS 1, 2 4, An 7/ Si 

= — (7/ + a/i +p)S-\ 

at 7 7 



A/i + /i + a 2 (/! 3 + 3/17), (79) 
2A7 + 2h 1 w-fS + 12a V7 + 7i + «V + P 2 , (80) 



+ (^f)(l + ZS)(l-S), (81) 

Equations (179 i) -( l8Ti) clearly show that fi(t) is coupled with j(t) and S(t). Figure 
10(a), 10(b) and 10(c) show time courses of /i(t), j(t) and S(t) given by Eqs. (179]) - 
(18T1) when mean-driven [Eq. (l60ll with A = 0.4, A), = 0.1], fluctuation-driven [Eq. 
( 16"T1) with B = 0.2, Bf, = 0.05] and synchrony-driven pulse inputs [Eq. (1621) with 
C = 0.8, Cb = 0.2], respectively, are applied to a neuron ensemble with a = 0.35, 
(3 = 0.1, u> = 0.5 and iV = 100. We note in Fig. 10(b) that /i(t) is slightly increased 
when 7/(t) is increased at 40 < t < 60, while fi(t) is independent of 7/(i) for 
F(r) = — Xr and G(r) = r [Fig. 6(b)]. A similar increase in /i(t) is observed when 
synchrony-driven input of Si(t) is applied as shown in Fig. 10(c), while no changes 
in Fig. 6(c). A peculiar time dependence is obtained in 7(t) for the mean-driven 
input in Fig. 10(a) compared to that in Fig. 6(a). Except these points discussed 
above, responses in the case of G(r) = r 2 shown in Fig. 10 are similar to those in 
the case of G(r) = r shown in Fig. 6. 

A number of neuronal experiments have not reported a systematic changes in 
firing rates while the synchronization within an area is modulated [6] . In particular, 
the synchrony is modified without a change in firing rate in some experiments [7J, 
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[9j [TT]. It has been pointed out that such phenomenon may be accounted for by a 
mechanism of a rapid activation of a few selected interneurons [12]. Recently the 
absence of a change in firing rate is shown to be explained if the ratio of excitatory 
to inhibitory synaptic weights of long-range couplings is kept constant in neuron 
ensembles described by the IF model [33] • It would be interesting to examine various 
cases of F(r) and G(r), changing their functional forms, which is left for our future 
study 

5 Conclusion 

We have studied responses of neuronal ensembles to three kinds of inputs: mean-, 
fluctuation- and synchrony- driven inputs, applying the AMM and DSs to the gener- 
alized rate-code model [34J . The ICA analysis of our results has suggested that mean 
rate, fluctuation (or variability) and synchrony may carry independent information. 
It would be interesting to examine this possibility by neuronal experiments using in 
vivo or in vitro neuron ensembles. 

One of advantages of our rate-code model given by Eqs. ([I])-© is that we can 
easily discuss various properties of neuronal ensembles, by changing F(r), G(r) and 
H(u), and model parameters. We hope that our rate-code model shares advantages 
with phenomenological neuronal models such as the Hopfield [60J and Wilson-Cowan 
models [61J. Although our rate-code model has no biological basis, it might be 
derived from spiking-neuron models with approaches previously adopted in Refs. 
[56j-[59j. It is promising to apply the generalized rate-code model to more realistic 
neuronal ensembles with excitatory and inhibitory dynamical synapses. 
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Appendix A: Derivation of AMM equations given by Eqs. 

(HSD-dTHD 

Using Eqs. (□), (EJ) and ©, we get 

(IT ■ 

= F( ri ) + H( Ul )+5I t (t)+ 0(^)^(1)+ Ut), (Al) 

with 




The Fokker-Planck equation for the distribution of p{{rk},t) in the Stratonovich 
representation is given by [H21 E3] 

|p(K},t) = -Y.-^{{ F (r i )+H(u t )}p({r k },t)} 
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+ \ E E ^rVSni + ^ + 0(i - <%)] P({r fc }, t)L 

+ ^T.^{G{r l )^-[G{r i ) P {{r k }^t)]}. (A3) 
Equations of motion for moments, (r») and (r$ r^), are derived with the use of FPE 

' ,{r,) (F(r i ) + if( Mj )) + ^(G / (r i )G(r l )), (A4) 



dt 



dt 2 

= (r, [Ffa) + #(«,)]) + fa [F(r 4 ) + H{ Ul )}) 



a 2 



+ y [(r ? G"(r J )G(r,)) + (r J G"(r l )G(r J ))] 



+ %[a 2 (G(r i ) 2 ) + 7/ + /3 2 ] + (l-^-)C/- (A5) 
We may obtain Eqs. (fT6 i) -( fl8{) . by using the expansion given by 

^ = fi + Sri, (A6) 



and the relation given by 



dfj. 1 x d(rj) 

dt ~ iV^^T' ( j 

d 1 l dJXSnft / a n\ 

dt = ivE^^' ( ) 

dp 1 d(5r^r^ 

For example, Eq. ffl6l) for dfi/dt is obtained as follows. 

^E(^W) = /o + M, (A10) 

^EW) = (ah) 
^E^VOGN) = ^1 + 3(^3 + ^2)7- (A12) 

Equations (fTTl) and ffl8|) are obtainable in a similar way. 

Appendix B: Jacobian matrix of Eqs. (1471)- (j49l) 

It is better to adopt the basis of (p, 7, p) than to adopt that of (p, 7,6'), in 
making a linear stability analysis. In the former basis, equations of motion given by 
Eqs. (TIfil) - (TTBl) become for F(r) = — Xr and G(r) = r: 
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d>y , 2h x wN 



± = -2\p + 2h 1 wp + ^(l + ZS I ) + 2a 2 p+^(a 2 p 2 + {3 2 ) 



(B2) 
(B3) 



The stability of the stationary state may be examined by the Jacobian matrix 
of Eqs. (B1)-(B3). With the use of cu = C13 = C32 = in the matrix, we get its 
eigenvalues given by 

a 2 

Ai = c n = -A + — + h lW , (B4) 

2h x w 



A 2 
A 3 



C22 
C33 



-2A + 2a z - 



-2A + 2a 2 + 2h x w, 

from which the stability condition given by Eqs. ( 1581) and ( 1591) is obtained. 



(B5) 
(B6) 



Appendix C: Quantities in temporal and rate codes 

We will discuss the relation between quantities calculated in the temporal code 
and those obtained in the rate code. In the former, the interspike interval (ISI) 
defined by 

T?(t) 



±k 1 k 

''m+l ^in> 



for t\ n E [t- At/2,t + At/2), 



(CI) 



plays an important role, where t\ n denotes the nth firing time of neuron % for a given 
fcth trial. When we assume that the firing rate is given by the inverse of ISI averaged 
over M trials: 



M 



(C2) 



k=l 



the ISI histogram (ISIH) is given by 



dr 






' dT 







(C3) 

where p(r) denotes the distribution for firing rates. For example, in the case of 
F(r) = — Ar, G(r) = r and (3 = 0, we get stationary distributions given by [31] 

2H 



p(r) oc r-( 2A / Q exp 
P{T) oc T^'-^exp 



a 2 r 



(2H 



T 



(C4) 
(C5) 



where P(T) is expressed by the gamma distribution. 

The synchrony of ISI may be evaluated by its equal-time, auto and mutual 
correlations averaged over N neurons and M trials, as given by [6] [50]- [52] 



A T (t) 
M T (t) 



N 



iv 1=1 

1 



(C6) 
(C7) 
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with 

^)4E% (C8) 

i 

A simple calculation using Eqs. <\G2\i . (106|) and <\G7\i with the condition: 

| n-ii{t) (C9) 

M*) = ^ !>(*). ( C1 °) 



leads to 



> T M - 4& = ^r. (on) 



/x(*) 4 

where A(t) and M(t) denote auto and mutual correlations, respectively, for firing 
rates given by Eqs. (1391) and (I4UI) . After Eq. (|44p . we may define the normalized, 
temporal synchrony for ISI data, S T (t), given by 

S T (t) = (C13) 
Substituting Eqs. flCTH and flClH to Eq. (lCl3L we get 

S T (t) ~ S(t) = (C14) 

where S(t) stands for the normalized, rate synchronization for firing rates given by 
Eq. fl20D. 

In order to assess the synchrony, one calculates the variance of ISI given by 

CUt) = (C15) 

From Eqs. (jC2| . flC9|) . flCTol) and (jCI5|) . we get [M] 



C#(f) ~ CV(f) = (C17) 

where Cy(t) is the rate variance given by Eq. (|24p . Stationary distributions given 
by Eqs. (JS1D and yield 

Cy = . ° (C18) 
>/2(A - « 2 ) 

^ = ° (C19) 
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which satisfy Eq. (1C17D for a 2 < A. 

According to neuronal experiments, spike train variability is correlated with lo- 
cation in the processing hierarchy [SI]. A large variability has been observed in 
hippocampus (Cy ~ 3) [SS] and in visual VI and MT of monkeys (Cy = 0.5 ~ 1.0) 
[66J whereas a variation of motor neurons is very small (Cy = 0.05 ~ 0.1) [ST]. For 
an explanation of the observed large Cy, several hypotheses have been proposed: 
(1) a balance between excitatory and inhibitory inputs [68j[69j, (2) correlated fluc- 
tuations in recurrent networks [70J, (3) the active dendrite conductance [71] . and (4) 
a slowly decreasing tail of input ISI of T~ d (d > 0) at large T [72] . Equations ( 1C17I) - 
(!CT9l) show that the variance is increased with increasing multiplicative noise, which 
may be an alternative origin (or one of origins) of the observed large variability. 
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Figure 1: (Color online) Stationary values of fi Cy and S as a function of [ij for 
w = 0.0 (dashed curves) and w = 0.5 (solid curves) with 7/ = 0.2 and Sj = 0.1 
(A = 1.0, a = 0.0, (3 = 0.1, and iV = 100), the chain curve expressing \i = jij. 



Figure 2: (Color online) Stationary values of Cy and S as a function of Cyj for 
w = 0.0 (dashed curves) and w = 0.5 (solid curves) with fij = 0.2 and Sj = 0.2 
(A = 1.0, a = 0.0, (3 = 0.1, and N = 100), the chain curve expressing Cy = Cyi- 



Figure 3: (Color online) Stationary values of \i and S as a function of Si for w = 0.0 
(dashed curves) and w = 0.5 (solid curves) with fij = 0.2 and 7/ = 0.2 (A = 1.0, 
a — 0.1, (3 — 0.0, and iV = 100), the chain curve expressing S = Si. 



Figure 4: /1 as a function of \ii for various a values with A = 1.0 and w = 0.0. 



Figure 5: S as a function of Si for various iV with /i/ = 0.1, 7/ = 0.0, A = 1.0, 
a = 0.5, (3 = 0.2 and w = 0.5. 



Figure 6: (Color online) (a) The time courses of //(£), j{t), S(t) and Cy(t) for 
mean-driven pulse input given by Eq. (|60!) with A = 0.4, A 6 = 0.1, 7/ = 0.1 
and Si = 0.1: (b) those for fluctuation-driven pulse input given by Eq. (loT|) with 
-B = 0.2, B b = 0.05, Hi = 0.1 and S7 = 0.1: (c) those for synchrony-driven pulse 
input given by Eq. with C = 0.4, C b = 0.1, /i 7 = 0.1 and 7/ = 0.1. Solid 

and dotted curves express results of the AMM and DS, respectively: chain curves 
express inputs of /i/(t), 7i(i) and Si(t) (A = 1.0, a = 0.1, (3 = 0.1, w = 0.5 and 
N = 100). 



Figure 7: (Color online) (a) The time courses of /x(t), j(t), S(t) and Cy(t) for 
mean-driven sinusoidal input given by Eq. (15^]) with 7/ = 0.1 and Si = 0.1: (b) 
those for fluctuation-driven sinusoidal input given by Eq. (1641) with hi = 0.1 and 
Si = 0.0: (c) those for synchrony-driven sinusoidal input given by Eq. fl65|) with 
fii = 0.1 and 7/ = 0.1. Solid and dashed curves express results of the AMM and DS, 
respectively: chain curves express inputs of Hi(t), ji(t) and Si(t) (A = 1.0, a = 0.1, 
= 0.1, w = 0.5 and N = 100). 
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Figure 8: ICA separation of the AMM result for mean- and synchrony- driven inputs: 
(a) source input signals: (b) output signals: (c) extracted signals by fast ICA [55] 
(A = 1.0, a = 0.1, /3 = 0.1, w = 0.5 and N = 100)(see text). 



Figure 9: ICA separation of the AMM result for mean-, fluctuation and synchrony- 
driven inputs: (a) source input signals; (b) output signals: (c) extracted signals by 
fast ICA [53] (A = 1.0, a = 0.1, (3 = 0.1, w = 0.5 and N = 100) (see text). 



Figure 10: (Color online) (a) The time courses of /i(t), 7(t), S(t) and Cyif) for 
mean-driven pulse input given by Eq. (I6"0l) with A = 0.4, A b = 0.1, 7/ = 0.2 
and Si = 0.2: (b) those for fluctuation-driven pulse input given by Eq. (loT|) with 
B = 0.2, B b = 0.05, \ii = 0.2 and Sj = 0.2: (c) those for synchrony-driven pulse 
input given by Eq. with C = 0.8, C h = 0.2, ^ = 0.3 and 77 = 0.2. Solid and 
chain curves express results of the AMM and inputs signals, respectively, in the case 
of F(r) = -Ar and G(r) = r 2 (A = 1.0, a = 0.4, (3 = 0.1, w = 0.5 and N = 100) 
(see text). 
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